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Abstract. We start our lectures with introducing and discussing the general notion 
of multifractality spectrum for random measures on lattices, and how it can be 
probed using moments of that measure. Then we show that the Boltzmann-Gibbs 
probability distributions generated by logarithmically correlated random potentials 
provide a simple yet nontrivial example of disorder-induced multifractal measures. 
The typical values of the multifractality exponents can be extracted from calculating 
the free energy of the associated Statistical Mechanics problem. To succeed in such a 
calculation we introduce and discuss in some detail two analytically tractable models 
for logarithmically correlated potentials. The first model uses a special definition of 
distances between points in space and is based on the idea of multiplicative cascades 
which originated in theory of turbulent motion. It is essentially equivalent to statistical 
mechanics of directed polymers on disordered trees studied long ago by B. Derrida 
and H. Spohn in |12| . In this way we introduce the notion of the freezing transition 
which is identified with an abrupt change in the multifractality spectrum. Second 
model which allows for explicit analytical evaluation of the free energy is the infinite- 
dimensional version of the problem which can be solved by employing the replica trick. 
In particular, the latter version allows one to identify the freezing phenomenon with 
a mechanism of the replica symmetry breaking (RSB) and to elucidate its physical 
meaning. The corresponding 1-step RSB solution turns out to be marginally stable 
everywhere in the low-temperature phase. We finish with a short discussion of recent 
developments and extensions of models with logarithmic correlations, in particular in 
the context of extreme value statistics. The first appendix summarizes the standard 
elementary information about Gaussian integrals and related subjects, and introduces 
the notion of the Gaussian Free Field characterized by logarithmic correlations. Three 
other appendices provide the detailed exposition of a few technical details underlying 
the replica analysis of the model discussed in the lectures. 
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Investigations of multifractal measures of diverse origin is for several decades a very active 
field of research in various branches of applied mathematical sciences like chaos theory, 
geophysics, oceanology, climate studies, and finance, and in such areas of physics as turbulence 
and statistical mechanics [IJ, and theory of quantum disordered systems [2]. The main 
characteristics of multifractal patterns of data is to possess high variability over a wide range 
of space or time scales, associated with huge fiuctuations in intensity which can be visually 
detected (see fig. 1). Another common feature is presence of certain long-ranged powerlaw- 
type correlations in data values. 




Figure 1. Multifractal probability density for a model of quantum particle at 
the critical point of Quantum Hall Effect, see [2j. Courtesy of F. Evers, A. 
Mirlin and A. Mildenberger, unpublished. 

To set the notations, consider a certain (e.g. hypercubic) lattice of linear extent L in 
A^— dimensional space, with M ~ standing for the total number of sites in the lattice. 
The measures of interest are usually defined via weights pi associated with every lattice site 
i = 1, 2, . . . , M and appropriately normalized to the total weight equal to unity as sketched 
below: 




Figure 2. A square lattice with weights attributed to the lattice sites. 

One can imagine a few different spatial arrangements of weights pi across the lattice 
sites. In the case of simply extended measures the weights are of similar magnitude at each 
lattice site, the normalisation condition then implying the scaling pi ~ in the large- 

M limit. As a generalisation of the above example one can imagine the non-zero weights 
Pi supported evenly on a fractal subset of lattice sites of effective dimension < N^j < N . 
In the limiting case of N(,f = we then deal with localised measures characterized by the 



Multifractality and Freezing 



3 



weights Pi essentially different from zero only inside one or few blobs of finite total volume. In 
such a situation weights stay finite even when M — t- c«, that is pi = 0(M°). Finally, in the 
most interesting case of multifractal measures the weights scale differently at different sites: 
Pi ~ M~"' [jj] The full set of exponents < Oj < oo can be conveniently characterized by 
the density p{a) = Yli=i ^{o^ ~ (^i) whose scaling behaviour in the large-M limit is expected 
to be nontrivial: p{a) ~ M^^"\ with the convex function /(a) known in this context as 
the multifractality spectrum or singularity spectrum, see Fig. 3. In view of the identity 
p{a) da = M we see that at the point of maximum a = qq we must have /(ao) = 1- 
Note also that the total number m(a) = p{a) da of sites of the lattice characterized by 
the scaling exponents Oj < a(< oq) satisfies for M ^ 1 the inequality m{a) ~ M'^*^"-* > 1, 
hence f{a) > for a < a^. Modifying this argument one can show f{a) > also for a > oq- 
The condition f{a) = defines generically the minimal a_ and maximal a+ threshold values 
of the exponents which can be observed in a given typical pattern. Note that the constraint 
Pi < 1 implies a_ > 0. 




Figure 3. Shape of a typical multifractality spectrum. 

An alternative, frequently more practical way of describing multifractality is via the set of 
exponents Tg characterizing the large-M behaviour of the so-called inverse participation ratios 
(IPR's) Pq which are simply the moments of the corresponding measure: 

M „ 

1=1 

Substituting in the above definition the relation p{a) ~ M^^"^ one can evaluate the integral in 
the large-M limit by the the steepest descent (also known as Laplace) method, see Appendix 
A. One then finds the relation between Tq and f{a) given by the Legendre transform: 

Pq^M-^", Tq = qa-f{a) where 'Z = £ • (2) 

In particular, at the point of maximum q = and as from the very definition tq = — 1 we 
immediately see that /(ao) = maxQ{/(Q)} = 1, cf. Fig. 3. 

II Usually one defines exponents via the relation pi ^ x^-Nai j_g_ j^-^q reference to linear scale L 
instead of the total number of sites M ~ L^. We however find it more convenient to get rid of trivial 
spatial dimension factor N , and concentrate only on essential parameter behaviour. 
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The above description is valid for multifractal measures of any nature. In recent years 
important insights were obtained for disorder-generated multifractality, see [2] and [3] for a 
comprehensive discussion in the context of Anderson localisation transitions, and [H |5] for 
examples related to Statistical Mechanics in disordered media which are closer to the context 
of the present lectures. One of the specific features of multifractality in the presence of disorder 
is a possibility of existence of two different sets of exponents, Tq versus fg, governing the scaling 
behaviour of typical Pq versus disorder averaged IPR's, < Pq >~ M""^'. So by definition 

InM ' ^" InM ' 
Here and henceforth the brackets stand for the averaging over different realisations of the 
disorder. The first type of averaging featuring in the above equation is traditionally called 
in the literature "quenched" , and second one is known as "annealed". It is known that the 
"quenched" values correspond to values of exponents which one finds in a "typical" realisation 
of disorder. The possibility of "annealed" average to produce results different from typical 
is related to a possibility of disorder-averaged moments to be dominated by exponentially 
rare configurations in some parameter range. A related aspect of the problem is that the 
"annealed" multifractality spectrum recovered from the multifractal exponents fq via the 
Legendre transform ([T]) can be negative: f{a) < 0, see fig. 4. 




Figure 4. Shape of an "annealed" multifractality spectrum with negative parts 
(dotted) extracted from the disorder-averaged moments and reflecting exponentially 
rare events, see the text. 

Indeed, those values reflect events which are exponentially rare [6] and need exponentially 
many realisations of disorder to be observed experimentally or numerically. On the other hand, 
as was noted in [2|, when dealing with typical multifractality spectrum /*^^(a) by exploiting 
the relation ([1]) one has to specify the limits of integration over a to be precisely a_ < a < «+. 
IPR moments are then given by 

ptyp ^ n M-«"+-^*'"'(")(ia ~ M-^"' , (4) 

and calculating the above integral by the steepest descent method reveals that typical (that 
is quenched) exponents Tq = r^^^ are related to py^{ct) by Legendre transform only in the 
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range -^\a+ = Imin < 9 < (Imax = ^U-) whereas outside that interval the exponents behave 
linearly in g, that is Tq'^ = qa±, see fig. 5. We will not dwell on the differences "quenched" 
vs. "annealed" exponents further and direct the interested reader to the recent works [5] and 
[3] for more detail and further references 0. In the present set of lectures we will concentrate 
exclusively on calculating typical (=" quenched") values of IPR exponents for some class of 
models. 





r — ^typ 

1 ~ 'q 

— yf^'a — "'^ 


Qmin 

-L 


/ ^ Hmax 



Figure 5. g-dependence of typical ("quenched") multifractality exponents Tq. Dotted 
lines show linear behaviour, see the text. 

Introduced through the moments involving summation over all the lattice sites, cf. ([I]), 
the multifractality by itself says nothing about more delicate questions, for example about 
spatial correlations between weights at two different sites of the lattice with coordinates, say, 
xi and X2, separated by a given distance |xi — X2|. The most natural assumption which 
is satisfied by vast majority of multifractal measures of actual experimental interest is the 
power-law decay of correlations implied by full statistical spatial self- similarity of the random 
measure: 

(/(xi)p^(x2)) oc L-^2'(9,s)5-iV^(g,s)^ 5 = |xi - X2| . (5) 

As statistical homogeneity of the random measure implies for local averages {p'^{xi)) = 
'M Si^i Pi ~ L~^~^'^i the equation ([5]) after setting 5 ^ L yields the relation for exponents: 

s) + s) - 2 = + (6) 

which follows from assuming the decoupling (p'?(xi)p*(x2)) ~ (p''(xi)) (p*(x2)) at large 
separations 5 = |xi — X2I ~ L — )• 00. On the other hand, for sites separated by a single 
lattice spacing (5 = 1 we must have (p''(xi)p''(x2)) ~ (p'^''~^(xi)) ~ L-^-^'^q+a ^ which after 
comparing with ^ allows one to relate the exponents governing the spatial correlations to 
the multifractality exponents as [7] 

y{q,s) = l+ Tq+s, Z{q, S) = I + Tq + Ts - Tq+s (7) 

^ Note that unfortunately the definitions of the termination of the multifractality spectrum used in [5] 
and in [3] are essentially different. The work [S] uses the definitions set up in the comprehensive review 
[2] which could be consulted in case of confusion. 
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Further, it turns out to be instructive to exploit (0) for evaluating the following correlation 
function: 

(lnp(xi)lnp(x2)) = ^ (p''(xi)p''(x2)) |,=.=o (8) 

Remembering tq = —1, and (Inp(xi)) = — A^lnL^|g=o we obtain after straightforward 
manipulations the following fundamental relation 

(Inp(xi) lnp(x2)) - (Inp(xi)) (lnp(x2)) = -g^ In ^i^, ^2 ^ ^^!|±£| ^^^^ > q (g) 

valid for arbitrary self-similar multifractal field. In other words, we have demonstrated that 
multifractality plus statistical selfsimilarity and homogeneity of the random weights necessarily 
imply that the logarithms of such weights must be correlated logarithmically in space. 

Inverting such an argument suggests that possibly the simplest way to generate random 
multifractal weights in the lattice is by constructing quantities Inpj at every lattice site i 
as Gaussian-distributed random variables correlated in precisely the way prescribed by ([9]). 
The resulting model has a very natural interpretation in terms of the equilibrium statistical 
mechanics. Indeed, consider a single classical particle subject to a random Gaussian potential 
y(x). It is the standard fact of theory of random processes[8] that if such a particle moves 
under the influence of the thermal white noise according to the Langevin equation 

d 



i = -^V (x) + e(x, t), ^xi, tl)?(x2, t2) = 2T5{ti - t2) 

then the probability ^(x, t) to find such a particle at a point x of the sample of finite size L 
will converge to the equilibrium Gibbs-Boltzmann measure 

P(x, t oo) pf}{yi) = exp -/3V{x) 

characterized by the inverse temperature /3 = 1/T. The normalization Jj^|^^p^(x)(ix = 1 
implies the value of the partition function to be given by 



Z{(3) = / exp -(3V{x) dx . (10) 

J|x|<L 

As obviously Inp^(x) = const — l3V{x) the weights p^(x) according to our discussion will 
be multifractal if the potential V^(x) is chosen logarithmically correlated in space: 



(y(xi) y(x2)) = -5'in 



(xi - X2)^ + 
L2 



a<L, xGM^, (11) 



where we assumed |x| < L, and the parameter a stands for a small-scale cutoff. 

According to the general discussion, the multifractal structure of the Gibbs-Boltzmann 
measure can be extracted from the knowledge of moments 

P,= [ p^(x)dx = -|^~L-^-^ as L^oo. (12) 

Identifying M ~ {L/a)^ , the Eqs.([l2D and ([10]) imply the following expression for the typical 
exponents Tq in terms of the appropriately normalized free energy of the system 

Tq = \q\mqW)-Qmf^), -^(/5) = - Jirn^ ii^l^M . (13) 

As shown in the Appendix A, the most natural random field with logarithmic correlations 
corresponds to the so-called Gaussian Free Field (GFF) in two spatial dimensions N = 2, as 
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well as its one-dimensional subsets. It is one of the fundamental objects in physics and various 
issues of its statistics attracted a lot of interest recently in conformal field theory, Schramm- 
Loewner evolution, and two-dimensional quantum gravity, see e.g. some discussion in [9|. 
Technically the problem of extracting the multifractality exponents Tq for the GFF amounts 
to ability to calculate efficiently the disorder average of the free energy (fT3|) . Such task is 
in general considered to be one of the most difficult problems in the statistical mechanics of 
systems with quenched disorder and we will not be able to perform such calculation explicitly 
in A'^ = 2 GFF case 0- Instead, we are going to outline such calculation for two particular 
choices of the models with logarithmically correlated potentials where such calculation is indeed 
feasible. The first model uses a special definition of distances between points in space and is 
based on the idea of multiplicative cascades which originated in the theory of turbulence, see 
e.g. discussion and further references in In fact, the model is essentially equivalent to 
statistical mechanics of directed polymers on disordered trees studied long ago in the seminal 
paper by B. Derrida and H. Spohn [T2j. Our second model will use standard Euclidean 
distances but exploits high dimensionality of the embedding space: — )• oo. Although 
the details of the two models and the corresponding methods of solution may look rather 
different, there is a general consensus that they address essentially the same physics: the 
so-called freezing transition common to all disordered systems with logarithmic correlations. 
And indeed we shall see that the resulting multifractality spectrum will be identical. In the 
final section we will give a short account of recent works on different aspects of logarithmically 
correlated potentials. 

2. Statistical mechanics for logarithmically correlated potentials generated 
by multiplicative cascades 

The construction we are going to describe below can be easily carried out in any spatial 
dimension, but for simplicity we consider the one-dimensional case of an interval of length L 
with the left end at the origin. With each point < X < L of such an interval we can associate 
an infinite binary string generated by expansion 

X = l(^ + || + ... + |^ + ...) =(X1X2X3...X„...) (14) 

where each Xn is either or 1. For some numbers the binary string is not unique but by choosing 
the expansion with infinite number of zeroes to the right it can always be made unique ( e.g. 
we use for L/2 the string (100 . . .) rather than (0111 ...)). Then for any two points X and Y 
in the interval we can introduce the distance function defined as (i(X, Y) = where n is the 
maximal number of first binary digits shared by X and Y. For example, if X = (0 *****.. .) 
and Y = (l*****...) then n = 0, hence d(X, Y) = (which is obviously the maximal possible 
distance between the points in the interval), if X = (00 ****...) and Y = (01 ****...) then 
n = 1, hence (i(X, Y) = etc. One can check that such a function ci(X, Y) indeed satisfies 
all the axioms for the distances: (i) (i(X, Y) > 0,VX ^ Y, and d(X, Y) = implies X = Y 
(n) d(X, Y) = d(Y,X) and the triangle inequality (ifi) (i(X, Y) + (i(Y,Z) > (i(X, Z) for any 
triple X,Y,Z. 

+ Actually, in recent years some sophisticated probabilistic methods were developed which allowed to 
address somewhat similar questions for GFF, see e.g. [T^ and the references therein. That development 
however goes beyond the remit of the present lectures. 
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Now, let us associate with every point X an infinite set of random i.i.d. variables 
(pk(K), k = 0,1,2, ... ,00 with zero mean and variances chosen to satisfy: 

(X) (Y)) = 2/ In 2 Si,kS{xix2X3...xk),{yiy2y3-yk) (1^) 

where we used the Kronecker symbol: Sa,b = 1 for A = B and zero otherwise, for any two 
objects A and B of arbitrary nature. Finally, with any point X of the interval we associate a 
random potential V (X) according to the rule 

00 

y(x) = </,o(x) + <^i(x) + ... = ^,^fc(x) . (16) 

fe=0 

This construction implies for any X 7^ Y: 

00 

{V (X) V (Y)) = {^l (X)) = 2^2 ln2 (n + 1) , (17) 

fc=0 

where we assumed that the two points X and Y share precisely n first digits in the binary 
expansion. This implies that they are separated by the distance (i(X,Y) = hence the 

above formula takes the form 

(F(X)F(Y)) = -2<72ln^&^, X^Y. (18) 

Wc see then that with respect to the chosen distance the constructed random potential is 
logarithmically correlated in space. When dealing with logarithmically correlated potentials 
one has to ensure the proper rcgularization at small distances, as the logarithm obviously 
diverges for X ^ Y. Various regularization schemes are possible, and in the present situation 
one of the most natural is to replace continuous space of the interval with a discrete lattice 
structure. In the particular case under consideration we introduce a "lattice" of 2^ = M sites, 
each site located at one of the points X^v = = 0, 1, 2, . . . , 2^ — 1. We can visualise 

this construction via the tree diagram, associating the random fields (X) to every branch of 
the tree as sketched in Fig. 6 for i^T = 3: 





K = 3 






100 




<t>l{{) 






■#'2(00) / ^ 


02(01) ?>2(l(t/ \ 


02(11) 


<#'3{ooo^\ ... y 




\'#>3(111) 
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3 4 5 6 
8 8 8 8 


7 

8 



Figure 6. Lattice of 8 sites and the corresponding tree diagram associating random 
fields to every branch of the tree. 
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Now we can define the distances in the same fashion as before, but since the maximal 
number of common digits can be at most K we get for the variance of the random potential 
a finite value (cf. (jlip ): 

(X)> = In 2{K + 1) = -2^2 In | , (19) 

where we have introduced the effective lattice cutoff given by a = . For this regularized 

lattice version we can now introduce the well-defined Boltzmann-Gibbs weights 

P/3(X^) = -l-exp-/3y(X;v), ^i^(/3)= 5] exp-/3y(X;v) (20) 

and try to calculate the associated free energy (lnZ/f(/3)), hence to extract the multifractality 
exponents Tq, see (fT3]) . The value of the potential V{^n) associated with each lattice site Xat is 
obviously obtained by adding all the random fields (/'(X) along the unique path connecting the 
site to the top level of the tree diagram. This implies the essentially multiplicative nature 
of the cascade model for the weight factors exp{— /3y(X7v)}- The most efficient way to 
organize calculations amounts to exploiting such a multiplicative structure combined with the 
hierarchical organization of the model which is obvious from the tree diagram decomposition 
as shown in Fig. 7 below. The described structure implies that 



= e-'^*" [Z^kUP) + (21) 

where z']^^ {f3) corresponds to the left /right-hand subtree of the tree in Fig. 7 which is of the 
depth K — \ as, reflected in the lower index. 





1 00 


Left 


Right 





Figure 7. The tree diagram decomposition leading to recursive relations for the 
partition function. 

Note that the fields 0(X) entering zj^2i(/3) are statistically independent of those entering 

Zj^_i{(3). To make the direct use of the structure of the equation ([2T]) it is expedient to 
introduce the generating function 

G/^(p) = (e-^^^(^)), p>0, (22) 

which is simply the Laplace transform of the probability density of the partition function. 
Denoting the probability density of the distribution for the variable (po with Vi^cpo) and 
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exploiting that the variables Z^^\^{I3) and Z^^}_^{fi) are independent of each other and 
identically distributed (i.i.d.) we arrive at the relation: 

Gk{p) = I V{M (e-^-^^°[4"^.(/^)+4^^.('')l) #0 = / n^)Gl_, (pe-^*) #. (23) 



Precisely in the same way we can relate Gk-i to Gk-2, etc. in a kind of recursive 
procedure which starts with the obvious initial condition Go{p) = e~^. Finally, it turns 
out that the subsequent analysis becomes more transparent if one introduces a new variable 
X = — -glnp G (— cxDjCxd). We arrive therefore at the recursion relations 

Giix) = j V{^)Gf_^{x + ^) dcp, / = 1,2,...K and GoW = e"^""" , (24) 

where we have replaced Gi [p = e~^^) — )• Gi{x), with some abuse of notations. 

Note: If from the very beginning we had considered a tree with arbitrary constant 
branching s > 1 instead of the binary tree with s = 2 the above recursion would be simply 
replaced by 

Giix) = j V{<P)GU + <P) dcP, Go{p) = e-^"'^ . (25) 

where {V^ (X)) = 2g'^{K + 1) In s is the variance of the underlying logarithmically correlated 
potential, cf. ()19p . and M = is the total number of points in the lattice. 

To understand better the nature of the solution of the above equations in the 
thermodynamic limit — t- cxd it is instructive to consider the following limiting case for 
the branching parameter: s = l + 5, (5^1. This implies scaling the variable (j) in such a way 
that < (/>2 2^2 In s 2g'^5. 

To be specific, one may just wish to use the Gaussian distribution P(</>) = ^ exp — ■ 
Then the right-hand side of ()25p takes the form 

j ^7=^-3^ G]+', {x + <P)dc^^ 1^ e-4 Gj^f (x + gyV^) 

which after straightforwardly expanding in powers of 6 reduces (j25|) to 

(x) In Gi^i (x) + g'^—^Gi^i (x) 



Gi{x) = Gi^^{x) + 5 



dx"^ 



+ 0{5') . (26) 



Thus in such an approximation the function Gi{x) experiences only small change in one step 
of iteration: Gi{x) — Gi-i{x) oc 6. Introducing to this end the variable t = 15 and consider it 
to be continuous in the interval t G [0, tmax = K5 In M] we can replace Gi{x) — G(x, t) and 
approximately write to the leading order Gi{x) — Gi-i{x) k. 5-^G{x, t). In this approximation 
the relation (|26p is replaced by a partial differential equation on the function G{x,t): 

^=9'0 + GlnG, G(x,0)=e--^\ (27) 
We also note that (i) by its very definition the function G{x, t) satisfies the following conditions: 

< G(x, t) < 1, G(x ^ -00, t) = 0, G(x ^ 00, t) = 1 (28) 

and (ii) the values G{x,t) = and G{x,t) = 1 solves the equation (j27p . All these observations 
are typical for the partial differential equations having the so-called travelling waves solutions 
of the form 

G{x,t) = W[x - m{t)], ■^m{t) = c{t) ^ ct when t ^ 00, c = const > (29) 
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Figure 8. Sketch of a typical front of the travehing wave solution. 



where the constant c plays the role of the asymptotic velocity of the front propagation, see 
Fig.8. 

Substituting such a form to ([28]) . denoting t = x — m{t) (so that e.g. ^ = — ^^"T'(*)) 
we see that the partial differential equation in the limit t — )■ oo implies an ordinary differential 
equation for W{t) which can be written as: 

g^W + cW + ^U{W) = {), where U{W) = ^(hiW-^ (30) 

where we have introduced the notations W = and W = ^75-- Obviously, we can interpret 
the latter equation as the Newtonian equation describing the motion of a classical particle of 
mass g"^ on the interval of the fictitious "coordinate" W G [0, 1] in fictitious "time" r subject to 
the dissipative force ("friction") cW plus the potential force generated by the potential U{W) 
sketched in Fig. 9: 





U{W) 






unstable 







equilibrium 








W 


stable equilibrium 



Figure 9. Sketch of the potential driving the motion of a fictitious overdamped 
Newtonian particle, see the text. 

By inspection we see that the position W = Q corresponds to the maximum of the 
potential, hence it is unstable equilibrium, and = 1 is the stable equilibrium (minimum of the 
potential). As by its physical meaning TV < 1 the motion of such particle must be overdamped, 
that is it should approach the stable equilibrium W = 1 m. a. monotonic way as r — )• 00 (i.e. 
the damping should be strong enough to avoid oscillations around the stable equilibrium which 
would bring W out of the physical interval.) To this end, let us consider in more detail the 
motion in the vicinity of the stable equilibrium by expanding: W = 1 — v, v{t) ^ 1 so that 
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to the linear order VFlnVF ~ —v and (|3Up is reduced to the Hnear second-order differential 
equation g^v + cv + v = Q whose general solution is given by 

v{t) = Ae^+^ + Be^'\ X± = -^{c± ^/ c'' - Ag^) (31) 

To have a non-oscillatory ("overdamped") asymptotic behaviour for r — >• 00 is only possible 
for c > 2g, so that 

Be '\ c>2g ^32) 

Ere 3 , c = 2g 

To determine the value of c it is natural to recall that according to the definition (j22p 
Gk{x) = (exp{— e~^^.Z'i<-(/3)}) so that naively expanding for x — )• 00 gives Gk{x — )• 00) k, 
1 — e~^^ {Zxifi)) -|- . . .. Using the Gaussian distribution of the random potential chosen for 
the present model one finds 

{ZKm = M(e-P^^^^) = Mexpj^ < ^^(x >| « Me^'^'^r^M ^ ^i+z^^g^ ^ (33) 
which implies 

G/^(x ^ 00) « 1 - e-^^^-- , with Tmax =x- ct^^^, t^ax = In M and c = c(/3) = ^ + /^^^ -(34) 

The above formula for the velocity c = c(/3) ensures the consistency between the asymptotic 
behaviour in (|34|) and in (|32|) as for such a choice holds the relation (3 = -^{c — \/ — 4(7^). 
Moreover, the choice is also compatible with the condition for overdamped motion as c — 2^ = 
{/3g — 1)^ > 0. However there is some subtlety in that formula which is most apparent if 
we follow the function c(/3) starting from the high-temperature regime /3 <C g~^. We see 
that the wavefront velocity c(/3) decreases with increasing f3 (decreasing temperature) down 
to the minimal value c(/3 = g^^) = 2g, and then for T < Tc = g starts increasing again, as 
schematically shown below. 



v{t — )• cxd) 




a 



25 




Temperature dependence of the front velocity. The dotted branch is unphysical and should be replaced with the constant value c — 2g. 

A rigorous mathematical analysis of the travelling wave equations by Bramson|13] 
revealed that such conclusion is however not quite correct. Namely, Bramson proved that 
for the initial conditions of the type (j27p the actual velocity of the travelling wave front is 
indeed given by c{/3) = + 1^9^ for P < /3c = g~^, but sticks to the minimal value Cmin = 2g 
everywhere in the low-temperature regime (3 > (3c- Such a picture implies, in particular the 
asymptotic form W{t — t- 00) ~ 1 — e"^''^, or equivalently the asymptotics 

G{x, t)^l- e-/^-(^-^-™*), T <Tc = g (35) 
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Figure 10. 



so that the profile of the function G{x,t) turns out to be temperature independent ("frozen") 
everywhere in the lower-temperature phase. Such behaviour certainly signals of a kind of strong 
non-analyticity, as e.g. it invalidates the expansion of the exponent in (exp{— e~''^Zi^(/3)}) 
which underlay our "naive" analysis. It is therefore appropriate to call such a drastic change of 
the behaviour a phase transition, which is known in the literature as the freezing transition. 

Qualitatively, the same picture holds generically for an arbitrary branching s > 1, that 
is for the solution of recursive equation (|25p . Namely, an extension of the above analysis |12] 
shows that in the thermodynamic limit — )• oo the solution takes the form 

Gk{x) = W[x-c,K], c, = ||j^g^^^^(^^^_,^,^^^^ ^^^^ (36) 

where /3c is the point at which the function from the upper line in (j36p has its minimum: 
■^Ci3\p=f3^ = 0. Such a knowledge allows one to calculate our main object of interest, the mean 
free energy — /3F(/3) = limx-s>oo {^t^Zk{(3)). To this end it is convenient to use the following 
integral representation for the logarithm: 

r \e-P - e-P^] ^ . (37) 
Jo P 

Remembering the definition of the generating function: Gk{p) = {^e~^^^^'^)) and Gq{p) = e~P 
and also the relation p = e^^ we after averaging of (|37p arrive at the important identity: 

/oo 
[Go{x)-GK{x)]dx. (38) 
-oo 

Inspecting the travelling wave form of the solution (j36p we observe that in the limit K ^ 
the difference Go{x) — Gk{x) (sketched in Fig. 10) is approximately equal to unity inside the 
interval x G [■^,Kc^], and is negligibly small outside: 

This immediately produces the simple result for the limiting free energy: 

- /3F(/3) = hm 1 (In Zk (/3) ) = /3 (39) 

with cp given by (j36p . Remembering M ~ and using the relation ()13p for the typical 
multifractality exponents, we find 
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In particular, for the earlier considered case of the Gaussian distribution V{(t)) = 
exp — -At we find 

-mp) = ^'^F ''<* = '/^. (41) 

\ 2|:, /3>A = l/9 

We see that the only control parameter for the model is 7 = /3^//3^. After a simple calculation 
using ()40p we recover the multifractality exponents for this case, which we are going to present 
only in the range q> 1: 

(g-l)(l-7g), 0<7<^ 



'"(j>i 



9(1-^/7)', ^<7<1 • (42) 



0, 7 > 1 



The phenomenon of vanishing of the exponents rg>i in the low-temperature phase 7 = p'^g^ > 
1 is one of the manifestations of freezing. It is qualitatively interpreted in terms of the 
Boltzmann measure being essentially localised on a few sites for low enough temperature 
or strong enough disorder. The typical multifractality spectrum corresponding to the above 
exponents is obtained according to the Legendre transform ([2]) which gives 

where the expression in the first line formally assumes the range of exponents a- = (1 — ^/7)^ < 
a < 1 + 7 = oq, whereas in second line < a < 2.^/7 = oq. The upper bound oq here 
corresponds to the point of maximum of /(a) and is related to the formal restriction g > 1 in 
()42p . In fact however it is not difficult to find Tq for any q and show that the expressions ()43p 
are valid in a wider range a E [a_, a_|_] where the boundary a+ is the largest root of f{a) = 0. 




a 



multifractality spectrum multifractality spectrum 

in the high-temperature phase below freezing temperature 



Exploiting the relation (j4]) for the typical multifractality spectrum one has to specify 
the limits of integration over a to be precisely a_ < a < a+. Substituting there ()43p 
and calculating the integral by the steepest descent method reproduces the values (|42p of 
the quenched exponents, that is t^^^ = Tq. Such a calculation confirms that the change of 



Multifractality and Freezing 



15 



behaviour of the exponent r^^^ to hnear in q for 7 > 1 /q^ is induced by the dominance of the 
boundary point «_ in the integration over a, in agreement with general discussion after 

Thinking in terms of the multifractahty spectrum it is also easy to see that the freezing 
phenomenon at 7 > 1 is related to a__ = 0, when the leftmost end of the curve /(a) hits the 
vertical axis precisely at zero level: /(O) = 0, see Fig. lib. 



3. Statistical mechanics for logarithmically correlated potentials in 
Euclidean spaces of high dimensionality 

As was discussed in the Introduction, we consider the Gibbs partition function of a classical 
particle confined to a spherical box of some finite radius L. We denote the corresponding 
domain as {Dl : |x| < L}. As before our main goal is to calculate the ensemble average of 
the free energy 

F = -^lnZp, Z = [ exp -^V{x) dx , (44) 

where /3 = 1/T stands for the inverse temperature and dx is the standard volume element in 
A^— dimensional Euclidean space. The average of the logarithm of the partition function is one 
of the central problems in the whole physics of disordered systems, and is usually performed 
with the help of the so-called replica trick, i.e. the formal identity 



(InZ^) = lim-ln(Z^>, = £-^^"=1 ^('^•») TT dx, . (45) 



Dl 



a=l 



The random Gaussian-distributed potential V{x) is characterized by zero mean and the 
covariance specified by the pair correlation function (jlip . Performing the averaging over the 
Gaussian disorder in Eq. (05]) according to the formula ()A.9p , we in the standard way arrive at 
the following expression: 



(Z^) = e7"i-f / e ^^"^'"^ L ^ J n d^a , (46) 

■^Dl a=l 

where we recall the definition of the main control parameter of the problem: 7 = To 
achieve further progress one has to suggest an efficient way of working with the resulting 
multidimensional non-Gaussian integral. To this end one may notice that the integrand in 
Eq. (|46p in fact possesses a high degree of invariance: it depends on A^— component vectors Xa 
only via n(n + l)/2 scalar products Qab = XgXb, a < b, and is therefore invariant with respect 
to an arbitrary simultaneous 0(A) rotation of all vectors x^. Moreover, our choice of the 
integration domain respects this invariance. To this end, introduce A x n rectangular matrix 
X = (xi, ...,x„) such that the A— component vector Xj forms i—th column of such a matrix. 
Then the matrix Q = X'^X is n x n positive definite, whose entries are precisely the scalar 
products Qab = XaXb, a < b. An efficient method of dealing with integrals possessing such 
type of invariance is based on the fundamental identity 

/ ... / I{X^X) dxi...dxn=CN,n [ liQ) [detQ]^^^ dQ, (47) 

i|xi|<L i|x„|<L 

^(iv-li^) 

where CN,n = j^n-i N-k\ aiid we assumed A > n + 1. The integration domain in the 
llfe=o^V^~J 

right-hand side is simply = {Q > 0, qaa ^ L'^,a = l,...n}, the volume element is 

dQ = YlaKb^lab- The above formula seem to appear originally in IlMJ but has not been much 
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in use before it was independently rediscovered in the context of theory of random matrices 
in [15]. In [16] it was exploited in the present context. Since the relation turns out to be quite 
useful in a few applications we present in the Appendix B its derivation taken from [15| with 
the purpose of making the notes self-contained. 

Applying such a transformation gives in our case: 



•yn In — 

' ft 



,ln 



12 



D 



(Q) 



[det QY 



dQ , (48) 



So far all our manipulations were exact for any spatial dimension, provided N > n + 1. For 
any finite < oo no further simplifications seem possible, any ways to proceed to analysis of 
()48p are presently unknown and yet to be found. 

The situation is better if we agree to consider the dimension A^ as one more control 
parameter and let it to be large: N ^ 1. After appropriate rescaling of the coupling constant 
g — )• g^/N (i.e. 7 — )• A^7) and also rescaling the integration variables Q — )• we can rewrite 
the exact expression for the averaged replicated partition function in the following form 



(Zfi =CN,n ( yj 



^Nyn^ In - 



(detg) 



-(n+l)/2 -N<5>n{Q) 



dQ 



Dq 



where 



^n(Q) = -^ln(detQ)+7 5^1n 



a<b 



Maa + Qbb) - gab + 1 



(49) 



(50) 



and A^ is assumed to satisfy the constraint N > n. The final integration domain Dq is: 



Dq = {Q > 0, qaa < = 2L'/a^ 



l,...n}. The form of the integrand in Eq.([39 



is precisely one required for the possibility of evaluating the replicated partition function in 
the limit A^ — t- 00 by the multidimensional Laplace (also known as the "steepest descent" or 
"saddle-point") method. The effective free energy relevant for extracting the multifractality 
is then calculated by replica trick as (see (jl3p and (j45p ) 



m/3) 



lim 

M^oo 



InM 



— lim -$„(Q) 

L-s>oo in L n^O n 



lim 



(51) 



where we have replaced InM ~ A^lnL, and the entries of the matrix Q should be chosen to 



satisfy the extremal conditions: 
n{n + l)/2 equations: 



dqab 



for a < b. This yields, in general, the system of 



and 



[Q 



-11 



:{qaa + qbb) - g'ab + 1 



7 



{qaa + qbb) - qab + 1 



0, a 



0, a / 6 



1,2,, 



(52) 



(53) 



One should also ensure that the solutions to these equations respects the constraint 
qaa < for all a = 1,. . . ,n imposed by the presence of the boundaries of the integration 
domain Dq, and also the fact of Q being positive definite. However, the above equations 
obviously imply 

[Q-']aa = -Il [Q-']ab^ Va = 1, 2, . . . , n . (54) 

The above condition ensures that the matrix must have at least one zero eigenvalue 
(which corresponds to the uniform eigenvector with all components equal) which is obviously 
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inconsistent with constraints on Q. We interpret such a failure as manifestation of the fact 
that the functional ^n{Q) cannot achieve its extremum inside the domain Q > 0,qaa < R^- 
This means that such an extremum should be looked for at the boundary of the domain: 
Qaa = R^i Va = 1, 2, . . . , n. In turn, it means that when searching for such an extremum we 
only vary <5n(Q) with respect the off-diagonal entries, and therefore only have to satisfy the 
equation (f53]l . 

Our procedure of investigating the equations (I52I53P in the replica limit n — )• will follow 
the standard pattern suggested by developments in spin glass theory [18]. We first seek for the 
so-called "replica symmetric" solution, and then investigate its stability depending on 7. When 
the replica symmetric solution is found inadequate, it should be replaced by the hierarchical 
("Parisi", or " ultrametric" ) ansatz for the matrix elements Qab, with various levels of replica 
symmetry breaking. 



3.1. Analysis of the model within the Replica Symmetric Ansatz. 

The Replica Symmetric Ansatz amounts to searching for a solution to (|52p .()53 p within subspace 
of n X n symmetric positive definite matrices Q such that qaa = Qd = R^ ^ for any a = 1, . . . n, 
and qa<b = Qo, subject to the constraints < qo < R^ to ensure positive definiteness. Inverting 
such a Q yields the matrix of the same structure, with the diagonal entries all equal and 
given by 

V. = + '^"(^ - ^) f55) 

and all off-diagonal entries given by 



^° {R2-q,){R2+q^(^n-l)) ^^^^ 



Pd-PO = (^'^) 

-R^ - qo 



Note, that 



In the replica limit n — )• the equations (j53j) and (|56|) give in this way the equation for qo: 

{R'^-qoV (R'-qo + i) ^ ' 

It is convenient to define the variable do = R^ — qo satisfying < do < R^ and reduce ([58]) to 
the simple quadratic equation (7 + l)(io — do{R^ — 1) — i?^ = 0. Choosing the solution with 
do > and remembering that we are actually interested in the large— L limit R? = 2L^/a^ ^ 1 
we find 

^0 = ^7T^[^' - 1 + V(^2-l)2 + 4i?2(l + 7)] « (59) 
2(1 + 7) 7 + 1 

Now we should calculate the value of the functional ^n{Q) for the replica symmetric solution. 
It is easy to show that detQ = {R? — qo)"'~'^[R^ + (n — l)qo\, so that in the limit n — )• we 
easily find from ()50p 



lim -$„(Q) = In do - J In (do + 1) ~ -(1 + 7) In L + 0(ln a) (60) 

where we again considered the limit L ^ a. This shows that the effective free energy ()5ip is 
given by 

PTifS) = lim lim -$„(Q) = 1+7=1+ /3^g^ . (61) 
L-5>oo In L n-s>0 n 
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This coincides precisely with the high-temperature {T > T^. = g, i.e. 7 < 1) result for the 
"cascade model" of the previous section, cf. (|4ip . which is valid before the freezing mechanism 
becomes operative. Our next goal is to understand how the freezing emerges and is maintained 
for T < Tc. 

3.2. Analysis within the Parisi scheme of the replica symmetry breaking. 

The standard way of revealing the breakdown of the replica-symmetric solution is to perform 
a stability analysis following the pattern of the famous de Almeida-Thouless paper [TO] in the 
theory of spin glasses, i.e. magnetic systems with random interactions. Such analysis can be 
straightforwardly done for the present type of system, see Appendix D of the present lectures, 
and shows that for a given value of R the replica symmetric solution becomes unstable for 
the temperatures T < Tc = 9^-^- Therefore at low temperatures stable solution will have to 
be one with a broken symmetry in the replica space. To derive the corresponding expression 
for the free energy of our model we will follow a particular heuristic scheme of the replica 
symmetry breaking proposed originally by Parisi in the theory of spin glasses, see e.g. |18j . 
or more recently [20] 0. To make the present set of lectures self-contained we describe in full 
detail the structure of the matrix Q, the ensuing Parisi function x{q) and the main steps of 
the derivation in Appendix C in full detaij^. Here we just sketch those objects schematically 
for the convenience of the reader: 
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12b Step- wise Parisi function 
for finite integer n 



PI Schematic hierarchical structure of 

°' the matrix Q in Parisi parametrisation 



We are actually interested in the replica limit n — t- 0. According to the Parisi prescription 

* In recent years the use of the scheme was justified by alternative rigorous mathematical procedures. 
For the model under consideration the corresponding equations were re-derived recently by a rigorous 
methods in [5T] without any recourse to the powerful but ill-defined replica trick, 
jj The Appendix is taken verbatim from [TB], but we used this opportunity to correct the important 
formula (jC.12p which appeared in [IB] in a distorted form. 
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explained in detail in the Appendix C, in such a limit x{q) becomes non- decreasing function of 
the variable q and the system can be fully described in terms of such an object. The function 
depends non-trivially on its argument in the interval qo < q < qk-, with qo > Q and qk < qd- 
Outside that interval the function stays constant: 

x{q < qo) = 0, and x{q > qk) = 1. (62) 

In general, the function x{q) also depends on the increasing sequence of k positive 
parameters rrii satisfying the following inequalities 

< mi < m2 < • • • < mk < m^+i = 1 • (63) 

If the number of levels of the Parisi hierarchy K tends to infinity we may think of the function 
x{q) as continuous in the interval qo ^ Q ^ Qk, with possible jumps at the end of the interval: 
q = qo and q = qk- 

As is shown in the Appendix C, in the replica limit the following identity must hold for 
any differentiable function g{q): 

lim -Tr [g{Q)] = g (qa - qk) + ^ q' ( ^(^) ^q) dq . (64) 
"-^on Jo \Jg ) 

In particular, for the first term entering the replica functional Eq. (j50p application of the 
rule Eq. ([M|) gives 

1 f'i^ 1 

lim - [Tr In {Q)\ = In (g^ - g^) + / dq . (65) 

n^o n Jo " x{q) dq 



The last term in Eq. (j50p is also easily dealt with in the Parisi scheme (see Appendix C), where 
it can be written as 

k 

7 

=0 

by using explicitly the derivative of the generalized function Eq. ()C.6p . Using integration 
by parts and taking into account the properties Eq. (j62p we finally arrive at the required free 
energy functional for the phase with broken replica symmetry 



' lini V(mz+i - mi) In {qd - Qi + I) = - J ln{qd-q + l)x'{q) dq, (66) 



lim -$„(Q) = 



Ik I 

\n{qd-qk)+ / .g, (~\^~ dq 

qd-qk + jq x{q)dq 



- I (in {qd -qk + l)+ x{q) dq) = <P{x{q)} (67) 

The functional (f){x{q)} should be now extremized with respect to the non-negative non- 
decreasing continuous function x{q), whereas as we know the variable qd must be fixed to 
its boundary value qd = R^- To this end we find it convenient to introduce two parameters 
dmin = - qk, dmax = R'^ - qo Satisfying < dmin < dmax < R^ and also to use 
t = R? — q as the new integration variable simultaneously replacing (with some abuse of 
notations) x{q = R? — t) ^ x{t). Such a renamed function x{t) is now non-increasing in the 
interval t £ [dmin, dmax], and satisfies x{t < dmin) = 1, x{t> dmax) = 0. 

As the result, the above functional assumes a somewhat simpler form: 

dt fdmax x(t)dt 

- 2^{x{t)} = In {dm^n) + / "T— x + 7 In idm^n + 1) + 7 / (68) 

^rfmzn "mm + L - X{t)dt Jdrarn 
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Varying the functional Eg . (1681) with respect to such a function x{t) gives after due 
manipulations with integrals the expression 

fdmax 1 di 

-25^{x{t)}= S{t)6x{t) = Q, S{t)=-ij—- ^, (69) 

* dmin + /i . x{t) dr 

Requiring the variation to vanish therefore amounts to the condition S{t) = 0, Vt G 
[dmin,dmax]- As this obviously implies ^S{t) = we can differentiate Eg. ([69]) once, and 
immediately get the equation 

/"* t + 1 1 

dmin + / x{t) dr = -—, yt £ [dmin, dmax] x{t) = — (70) 

What remains to be determined are the values for parameters dmin and dmax- To this 
end, we substitute the value t = dmin into the first of relations Eg. ([70]) which shows that 

dmin — 7 ■ C^^) 

V7 - 1 

Next, we use the condition S{dmax) = 0, which in view of (169]) and x{t) = for t G [dmax,R^] 
gives the relation 



1 /-^^ di 



7 



dmax ~\~ 1 



J 2 • 



dm^n + ^ ™^ x{t) dr 



Substituting here the expressions (j70l7ip yields after a simple algebra the 7— independent 
result: 

i?^ — 1 

"'max — ^ ) V ' j 

completing the solution. According to the general procedure the solution makes sense as long 
as dmin < dmax, and using = g/T it is easy to check that the condition can be rewritten 
as T <Tc = g ^iZ\ which defines the low-temperature phase of the model for finite ii, with 
the same Tc as follows from the stability analysis (Appendix D). In the thermodynamic limit 
lim/j_^oo Tc = g, that is 7c = 1 as expected. 

The value of the functional at the extremum can be easily calculated by substituting 
x{t) = 7~i/^ for t G [dmin, dmax] and x{t) = for t G [dmaxjR"^] into ([68]) and using ([7T]) and 
(j72]) . This gives after some algebra 

- = In + ^ - ^(y^ - 1)2 In (V^ - 1) + 1(7 - 2^) In ^7 (73) 

which finally implies in the thermodynamic limit L — )• 00 for the effective free energy the value 

f3T{(3) = lim -i- hm -^n{Q) = 2 ^7 = 2/3ff • (74) 
L-s>oo m L n-s-O n 

In particular it shows that the free energy value in the low-temperature phase is frozen i.e. 
given by the temperature-independent constant T{f3) = 2g. This fully corroborates the picture 
obtained in the framework of logarithmic cascades of the previous section, see (j4ip . 

Before finishing this section it makes sense to discuss in more detail the picture associated 
with the freezing transition which manifests itself via the spontaneous breakdown of replica 
symmetry. The general interpretation of the freezing below Tc is that the partition function 
becomes dominated by a finite number of sites where the random potential is particularly 
low, and where the particle ends up spending most of its time [22]. For a more quantitative 
description of the particle localization, useful in the following, it is natural to employ the 
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overlap function defined as tlie mean probability for two independent particles placed in the 
same random potential to end up at a given distance to each other. Denoting the scaled 
Euclidean distance (squared) between the two points in the sample as P, and employing 
the Boltzmann-Gibbs equilibrium measure pp{^) = -^py exp — /3y(x) the above probability in 
thermodynamic equilibrium should be given by 

7r{V) = ([ dxip^(xi) / dx2P^(x2)5(p- J|xi-X2n ) (75) 

\-'|xi|<L J\JC2\<L V ^ J/y 

where again 6 denotes the Dirac's 5-function. The disorder averaging in (|75p can be calculated 
following the same standard steps of the replica approach as the free energy itself ( see 
Appendix A of |17j). With the function 7r(T>) in hand we can ask, in particular what is 
the probability for the particle in logarithmically correlated potential to end up at P = O(a^), 
i.e. at a distance of order of the small cutoff scale. The answer turns out to be zero in the 
high-temperature phase T > Tc, confirming the particle delocalization over the sample. In 
contrast, in the low-temperature phase T < T^. the probability is finite: vr (0(a^)) = 1 — T/T^, 
since both particles can be trapped by one and the same, or nearby favorable, deep minima. 
At a formal level such a behaviour is directly related to the shape of the function x{t) which 
in our case turned out to be rather simple and consisting of three flat regions (see Fig. 13b): 

x{0 <t < dmin) = 1, x{dmin <t < dmax) = , x{dmax < t < R^) = . (76) 

This essentially means that from the very beginning we could restrict ourselves to the first 
non-trivial level A: = 1 of the Parisi hierarchical scheme, see Eq. ijC.ll IC.2P instead of assuming 
the most general Parisi scheme for Q at the outset of our procedure. Such a simplified form 
(see Fig. 13a) of Q below the transition is typical for the random energy models and is known 
in the literature as 1-step RSB scheme, see e.g. [20] . The equilibrium values of the parameters 
go, qi and mi = m found from directly extremizing the corresponding functional ^<I>„((5)|n^o 
(or equivalently from solving the equations (j53p ) are given by 

qo = — - — ,qi=R — — -, and m=— . (77) 

2 - 1 Vl 

These values fully agree with those of the function x{t) obtained from the general Parisi Ansatz. 
Finally, in the Appendix D we discuss stability of the 1-step RSB solution for the logarithmic 
potential, and find it is actually marginally stable everywhere in the low-temperature phase. 
The latter feature is usually associated with the infinite-step Parisi Ansatz, see e.g. [20]. This 
is another manifestation of the fact that the logarithmic case is very special and shares both 
features of the full-scale infinite and 1-step replica symmetry breaking. 
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4. Summary, Historical background and Recent Extensions 

In this set of lectures we have addressed in some detail the spatial structures of the Boltzmann- 
Gibbs measure describing a single particle that thermally equilibrated under a random 
potential with logarithmic correlations. We have been able to calculate the multifractality 
spectrum of the measure and revealed the associated freezing transition by analysing the 
ensemble-averaged free energy in two special cases by two complementary methods. The 
first model introduced logarithmic correlations via employing the hierarchical "multiplicative 
cascades" construction and associated definition of the distance function. This way allowed 
us to perform the analysis of the freezing transition in the framework of a certain travelling 
wave equation satisfied by an appropriately defined generating function of partition function 
moments. In the second case the spatial dimension N of the system was assumed to be large 
which helped to employ the replica trick combined with the steepest descent method, and to 
relate freezing to the phenomenon of spontaneous replica symmetry breaking. In both cases 
the resulting free energy appears to be given by essentially the same expression. 

As the present day understanding of freezing and related phenomena has already a history 
of almost thirty years, it is certainly useful to be aware of a broader context of the problem 
under consideration. To this end it is appropriate to mention that an extreme "toy model" 
case of the problem in hand is represented by the famous Random Energy Model (REM) by 
Derrida where the freezing phenomenon was discovered and investigated for the first time 
[23\ I24j . The REM in some loose sense can be looked at as a limiting "zero-dimensional" 
A'^ = case of the model we studied elsewhere in this set of lectures. It amounts essentially 
to replacing the logarithmically correlated random potential by a collection of M uncorrelated 
Gaussian variables with the variances chosen to be scaled with M in the same way as in the 
logarithmic case: < V-^ >= 2g^hiM. REM is simple enough to allow explicit calculation of 
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the free energy by direct counting of degrees of freedom, and the result essentially coincides 
with ()4ip . A very informative account of the REM problem can be found in the fifth chapter 



Understanding quantitatively the generic statistical-mechanical behaviour of disordered 
systems for finite N is notoriously difficult, and even the simplest cases like our single-particle 
model still present considerable challenges. To this end we first need to mention a general 
attempt of investigating such model for finite dimensions < oo in the thermodynamic limit 
L — )■ c« undertaken in an insightful paper by Carpentier and Le Doussal [25j. The approach 
of Carpentier and Le Doussal was based on applying a kind of real-space renormalisation 
group treatment to the free energy distribution. The authors concluded that for finite spatial 
dimensions the model with logarithmically correlated potetial is really distinguished among 
others of similar kind. Namely, if correlations of the random potential grow faster than 
logarithm with the distance, then in the thermodynamic limit the corresponding Boltzmann- 
Gibbs measure turn out to be always localised at any temperature T < oo. At the same time, 
if the correlations decay to zero for large separations (such potentials are natural to call "short- 
ranged" ) than the Boltzmann-Gibbs measure turns out to be always trivially extended at any 
positive temperature T > 0. And only for the marginal situation of logarithmic correlations the 
true REM- like freezing transition indeed happens at some finite T = Tc > 0, at any dimension 
A > 1. Indeed, for that case the renormalisation group yielded after some clever albeit not 
fully controlled approximations a kind of travelling wave equation for the generating function, 
akin to (j27p . Fortunately, the logarithmic growth is not at all an academic oddity. The paper 
of Carpentier and Le Doussal can be warmly recommended for describing the present model 
in a broad physical context and elucidating its relation to quite a few other interesting and 
important physical systems, as e.g. quantum Dirac particle in a random magnetic field [26], 
and directed polymers on trees with disorder [12]. The latter works played the fundamental 
role in advancing the understanding of the freezing transition. Our presentation in the Section 
2 is actually based on an adaptation of material from |26j and [12], with the pedagogic example 
of branching tending to unity inspired by [2^ . 

Another line of research which deserves mentioning was pursued recently in [T7j where 
it was revealed that the picture of potentials with short-ranged, long-ranged, and logarithmic 
correlations presented in [25] is still incomplete, and misses a rich class of possible behaviour 
that survives in the thermodynamic limit L — )• oo. Namely, given any increasing function 
^{y) for < y < 1, it was suggested to consider Gaussian random potentials whose two-point 
correlation functions (covariances) take the following scaling form 



which generalizes our Actually, the above expression gives back (jlip for the special case 
^{y) = 5^(y ~ !)• As shown in \IL7\ the potential with the covariance (1780 can be constructed 
by a superimposing several logarithmically correlated potentials of the type (jlip with different 
cutoff scales a^, and allow those cutoff scales to depend on the system size L in a power-law 
way: Oj ~ < Ui < 1 . 

The equilibrium statistical mechanics of such system in the limit A — ?■ oo and L — )■ oo 
turns out to be precisely equivalent to that of the celebrated Derrida's Generalized Random 
Energy Model (GREM) see [28] and references therein. Namely, the system experiences a kind 
of freezing transition at the critical temperature Tc = \J ^'(1). Below this temperature the 




of i20j. 
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equilibrium free energy turns out to be in the thermodynamic hmit L — )• oo 

-^(r) = rz.,(r) + M^^il^^ + 2y^'y^(iy, o<r<T,, (79) 

where the parameter i/^, is related to the temperature T via the equation = $'(z/*). For 
T > Tc the free energy is instead given by 

_^(r) = T+li<ltJffl. (80) 

Using the two-point probability defined in (j75p these expressions for the free energy can be 
given a clear interpretation as describing a continuous sequence of " freezing transitions" which 
start at Tc and continue at all lower temperatures, with freezing happening on smaller and 
smaller spatial scales with decreasing temperature [17j . This is related also to the nature of 
the replica symmetry breaking, which requires for its description the full infinite sequence 
X — )• 00 of hierarchy levels in the Parisi scheme of Appendix C. Such a rich picture results 
in a more complicated multifractality spectrum /(a) which in contrast to (03]) is in general 
non-parabolic. However, it is appropriate to mention that the Boltzmann-Gibbs probability 
measures generated by the random potentials described in (j78p are rather peculiar, as for any 
non- linear function ^{y) they do not satisfy the standard spatial self- similarity property ([5]). 
Instead, it is easy to check that the exponents s) and s) governing the spatial decay 
of correlations between weights in ([5]) will be non-trivial functions of the variable ^"^^^^^^^ 
rather than simple constants. In this way, the exponents governing the decay of correlations 
for two points separated by the distance, say, |xi — X2I ~ L^/^ will be different from those 
separated by, say, |xi — X2I ~ L^/^. Though such behaviour is certainly not prohibited by first 
principles, it remains to be seen whether random multifractal measures with such peculiar 
spatial structure could appear in interesting applications in physics or other sciences. 

Although our lectures were centered around the notion of the multifractality spectrum, 
there is a different, and in essence deeper aspect of the freezing transition which attracted 
considerable research interest recently: the issue of the extreme value statistics [291 ISSl |30l [31] . 
This goes beyond the calculation of the ensemble-averaged value for the free energy F = 
—T\nZ{l3), but aims to describe precise form of the fluctuations around that mean value. 
Technically it amounts to our ability to calculate the shape of the generating function 
G{p) defined in (|22p in much finer detail (note that in the context of calculating typical 
multifractality exponents actual form of that function appeared to large extent irrelevant). 
As liniT-^oF = miux ^(x) it is obvious that at low enough temperatures the free energy 
fluctuations are dominated by the distribution of the deepest minimum of the random potential 
in a given sample. Classifying possible types of extreme value statistics for strongly correlated 
random variables is an open problem in probability theory with many important applications 
in natural sciences and beyond, see [25] and the references therein. In particular, it was argued 
in [25] that logarithmically correlated potentials represent a new universality class for extreme 
value statistics, and recent works |30t [3T] on extremes of the two-dimensional Gaussian free 
field (see definition of this important object in Appendix Al below) along various curves 
further substantiated that claim. Another aspect of the problem which certainly deserves to 
be mentioned here are intriguing but so far poorly explored connections to two-dimensional 
quantum gravity models as noticed in |34j . discussed in |25|, and most recently in [31j. Some 
speculations about relevance of the REM-type models in the string theory context can be 
found in [35 1. 

Finally, let us mention that there exists a completely different source of interest in 
multifractal random processes &: measures with logarithmic correlations motivated by growing 
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applications in financial mathematics, see e.g. |32j . |33j for the background information and 
further references. Although the questions addressed there are formally rather different, one 
can recognize a common mathematical structure. It is therefore natural to expect a fruitful 
merger of the two lines of research in the nearest future. 
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Appendix A. Elementary facts about Gaussian integrals and processes, the 
steepest descent method, and the Gaussian free field 

The fundamental role in applications is played by the standard Gaussian integral 

r e-f?/^+f'y ^ = 1 ei, Re{a) > 0, V6 (A.l) 
J-oo V2vr V« 

Suppose now we are interested in finding the asymptotic behaviour for large values of a 
parameter of the following integral 

e~^P^y)<P{x)dy, iV>l (A.2) 



where F{y) and (piy) are some given infinitely differentiable functions. It is clear that 
if the function F{y) is monotonically increasing/decreasing in the interval y G [2/1,2/2]! 
then the integral will be dominated by the vicinity of the left/right end of the interval, 
and for getting the leading asymptotics it is therefore enough to expand F{y) around the 
corresponding point up to the linear term only. For example for F'{y) > 0,Vy G [a,b], we 
write F{y) F{yi) + F' {yi){y — yi) + . . . which gives 

yi NF'[yi) 

where we assumed that generically 0(2/1) 7^ (otherwise one has also to expand ipiy) around 
2/1, which will change the result slightly). 

Similarly, if the function F{y) has a single maximum in some point 2/0 inside the interval, 
then subdividing the integration domain into two subintervals y € [2/1,2/0] and y G [2/0,2/1] we 
can apply the above consideration to each of the new intervals. For example, if F{yi) < F{y2) 
we have the same asymptotics as above in (IA.3p . whereas for F[yi) > F{y2) we have 

S/2 1 

e-^^(^)0(2/) dy « -— — e-^^(f2)</.(2/2)[l + 0(1/A^)] (A.4) 
j/i NF'{y2) 

Finally, the most interesting case arises if F{y) has a single minimum in some point 2/0 £ 
[2/1,2/2], that is F'{yQ) = and F"(2/o) > 0. In such a case the integral will be obviously 
dominated by the vicinity of the point of minimum, around which we can therefore expand as 
F{y) ~ F{yi) + ^—^p^{y — yif' + . . .. Substituting this approximation back to the integral and 
again assuming that generically (/>(2/o) 7^ we find after application of (jA.ip with a = F"{yQ), 
the asymptotics 



yi 



y2 / 9 

e-^^(^)0(2/) dy « e-^^(^")0(2/o)[l + 0{l/N)] (A.5) 
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These formulae represent the essence of the steepest descent (a.k.a. the Laplace) method of 
asymptotic evaluations of integrals. 

All the formulae can be naturally extended to the multivariable case. The multivariable 
generalisation of the Gaussian integral is given by 

(27r)'^/2 Vd^ ^ ' ' 

where n x n matrix A is assumed to be real symmetric Aij = Aji, and positive 
definite, i.e. all its eigenvalues Aj are positive. Then the inverse matrix A~~^ is well-defined 
and the determinant det A = ]Xi=i '^i 7^ 0. In fact introducing the scalar product for 
two vectors as (y, x) = j/jXj the quadratic form in the exponential can be written as 

AijUiUj = (y, ^y). The matrix is positive definite iff (y, Ay) > 0, Vy 0- 

The analogue of (jA.SP has the form 



-^^(y^"-y-U{y,, ...yn)dy,... dy^ « J .^^"^^"^ , e"^^(^i"-^")<A(yi, , . . . 2/„)|™n(A.7) 

V A*" detdo-F min 



where we assumed that the function F{yi, , . . . y„) has a single minimum at some point, and 

dpidyj 



f:j2 

mill is the n X n Hessian matrix {62F)ij — Ov dv -^(^i' ' ■ ■ ■ ^^i) evaluated at the point of 



minimum of F. 

Let us clarify the probabilistic meaning of the integral ()A.6p . Suppose that n real variables 
vi, . . . ,Vn are Gaussian-distributed, that is their normalized joint probability density of the 
vector V = . . . , is given by . . . , f„) = e~^^^'^'^K/ f^^d with some positive definite 



(27r)" 

matrix Aij. Denoting the averaging over such a distribution with the angular brackets (...) 
we can rewrite ()A.6P for any given vector b = . . . , 6„) as 

where the identities for the mean value and the pair correlation functions (a.k.a. covariances) 
immediately follow after expanding in the Taylor series with respect to 6j. 



Appendix A.l. Gaussian random fields: 



"massive" vs 



"free". 



The last expression is the basis for discussing properties of random processes (i.e. random 
functions V{x) of a single real variable x) which are a particular case of random fields 
representing random functions V{'x) of A^— dimensional vector x = (xi, . . . , xat). The field 
y(x) is called Gaussian if for any choice of the number n = 1, 2, . . . , cx) of points xi , X2, . . . , x„ 
in the space the joint probability density V{vi^ . . . ,Vn) of n values of the field in those points, 
that is vi = y(xi), f2 = ^(^2), ■ ■ ■ ■,Vn = l^(x„) are given by a Gaussian distribution with some 
matrix Aij. Such random field is uniquely determined by the two-point correlation function 



(the covariance) (y(xi)y(x2)) 



exp 



/(xi,X2) in terms of which the analogue of (jA.SP reads 

1 

2 



6(x)y(x) d^x 



exp 



/(X1,X2)6(X1)6(X2) dVf^''x2 



(A.9) 



for any suitable function 6(x). If we define the scalar product of any two functions a(x) and 
6(x) in the standard way as (a, b) = J a(x)6(x) d^x, we see that the quadratic form in the 
exponential of the right-hand side is (a, Fh) , where the linear integral operator F is defined via 



ff In fact the domain of validity of the formula (|A.6p is broader, and allows the matrix A to have 
complex eigenvalues with positive real parts 
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the kernel /(xi, X2). If one then defines the inverse operator as A = F ^ , the joint probability 
density of the random field V{'x) can be symbolically written using the scalar product as 



V [y(x)] = 1 exp 



-\{V,AV) 



(A.IO) 



where N is the suitable normalisation constant. 

To illustrate the latter approach, we briefly describe the paradigmatic example of the 
massive Gaussian field in N dimensions which is of importance for us here, and is also central 
for the modern theory of phase transitions. The probability of a given configuration V{'x) of 
such field is given by (jA.lOp with the quadratic form defined by 

(y, AV) = j (mV2(x) + ^2 [V1/(x)]2) d^x = j F(x) [n? - A] F(x) d^x (A.ll) 

where the "mass" m and the "stiffness" k of the field are two parameters, V is the gradient 
operator and A is the Laplacian: A = X^ILi Second form follows from the first one 
after applying the integration by parts and assuming that the random field V^(x) vanishes at 
infinity. In such an example the role of the operator A is obviously played by the second- 
order differential operator A = m? — A. Such operators are called "local" as their action 
on any function involves only values of that function and its derivatives in the same point 
of the space. Knowing A explicitly allows one to find the two-point correlation function 
(y(x) y(y)) = /(x, y) as the kernel of the operator inverse to A , hence satisfying the 
differential equation 

[m^-K^A] /(x,y) = 5(x-y), (A.12) 

where the Laplacian is assumed to act on the first argument, and (5(x — y) = J e^^^^^'^^ (27r)'^ 
stands for the appropriate Dirac delta-function. By applying the Fourier transform to the 
equation immediately gives the two-point correlation function as 

giq{x-y) ^JVq 

{m? + K2q2) {2-Ky 

To calculate the above integral it is convenient to use the identity (m^ -|- K^q^)^^ 
Jg°° e"*^™ ^'^ dt and change the order of integration, which gives 



00 







(y(x)y(y)) = / e-'^'dt j e^'i(-y)-^*'i^^ (A.14) 



1 dt _ 1 m^/2-ii^;v/2-i(^|x-y|) 



where we have used (jA.SP with Aij — )• 2K'^t5ij,h — )• (x — y) to evaluate the Gaussian integral 
in the first line, and K^{z) is the so-called Macdonald function, see the formula 3.471.9 of [36j . 
In particular, for N = 2 and m — )• we have from the expansion 3.471.9 of [36] 



y| « -,(A,15) 

m 



2K/m 

We conclude that the limit of 2D massless Gaussian field provides us with a random field 
with logarithmic correlations. 

The massless Gaussian field is also known in the modern literature as the Gaussian 
Free Field (GFF) and considered to be an object of fundamental importance. It can be 
defined on any domain D of A^— dimensional space using the following construction. Consider 
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an eigenproblem for the Laplace operator —A acting on functions in D, and denote ej(x), 
j = 1,2, ... ,00 its eigenfunctions corresponding to the Dirichlet boundary conditions ( i.e. 
vanishing at the boundary (9D) and let \j > be the corresponding eigenvalues. Then the 
functions ej(x) = --^ej(x) form an orthonormal basis of the Hilbert space H with respect to 

the so-called Dirichlet scalar (or "inner") product 

(/, g)= I (V/ • Vg) d^^ = - [ if- Ag) d^x (A.16) 

for functions /(x) on D vanishing at the boundary dD. Introduce now a set Q, j = 1,2, ... ,00 
of standard Gaussian independent, identically distributed real variables with mean zero and 
unit variance each: {(^j) = 0, = 1. Then the GFF V{x.) on the domain x G D is defined 
as the formal sum 

00 

y(x) = ^Oe,(x), (A.17) 

i=i 

from which it immediately follows that it is a Gaussian field with the covariance given by 



00 



(F(xi)F(x2)) = pe,-(xi)e,(x2) = - (A-^) (xi,X2) (A.18) 

j=i J 

which is nothing else but the Green function G(xi, X2) of the Laplace operator on the domain 
D. Note however that mathematically V^(x) is rather subtle (e.g. the sum in (IATTTI) does not 
converge pointwise and fails in general to be an element of the Hilbert space H ) . Because of this 
and other subtleties an extra mathematical care is needed to define the object fully rigorously, 
see references in [9]. The physicists however work with such an object without further ado, 
and we finish this section by two simple but important examples. Li the first example we deal 
with the GFF on a one-dimensional domain, the interval D = [0,1]. The Laplacian in one 
dimension is simply A = —-^ and the eigenfunctions/eigenvalues of the Dirichlet problem are 
given by e„(x) = ^/2 sin mrx, Xn — vr^fi^ so that the GFF in this particular case is given by a 
random Fourier series V{x) = Yl'^=i sinnvrx (compare with the periodic 1/f noise in the 
end of this Appendix). The corresponding Green function can be easily found to be given by 
G{xi,X2) = xi{l — X2) for X2 > xi and G{xi,X2) = X2{1 — xi) for X2 < xi. One immediately 
recognizes that the one-dimensional version of the GFF for such a domain coincides with the 
version of the Brownian motion called Brownian bridge, which is conditioned to return to the 
origin after a given time. 

Our second example is much more relevant in the context of the present lectures and 
deals with GFF defined on the two-dimensional disk: T) = \z\ < L where we use the complex 
coordinate z = x + iy. The Green function for the Dirichlet problem on such a domain is well 
known and is given by G{zi, Z2) = — ^ In xfei]^- particular, for any two points 1 2:1,2! ^ L 
(i.e. well inside the disk) the Green function reduces to expression equivalent to the full-plane 
formula (jA.lSp which is the basis for models with logarithmic correlations. 

Using the full-plane logarithmic GFF it is easy to construct various one-dimensional 
Gaussian random processes with logarithmic correlations. In particular, sampling the values 
of such GFF along a circle of unit radius with coordinates z = e**, t G [0, 2tt) we get a Gaussian 
process with the covariance {V{ti)V{t2)) = — ^ln|e**^ — e**^|. Such a process can be shown 
to be equivalent to a random Fourier series of the form V{t) = X^^i ["UnC*"* + Wne"™*] , 
where Vn,Vn are independent, identically distributed complex Gaussian variables with mean 
zero and variance {vnVn) = 1 (compare with the earlier Brownian bridge example). As the 
mean-square value (the "spectral power") of the coefficient in front of a given Fourier harmonic 
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with index n in this case decays hke 1/n such signals are known in many appHcations as 1// 
noises. 

Appendix B. Proof of the identity ( 1471) 

We start with identically rewriting the left-hand side of (j47p as 

j I {X'^X) dX = lim j I{Q)J,{Q) dQ (B.l) 

where 

J,(Q) = J e-f^'^I^'^lj (g - X^X) dX (B.2) 

and 5{x) stands for the appropriate Dirac (5— distribution in the matrix space. As usual 
(5— function can be expressed via the Fourier transform 6{x) = f e~*-^^^ its matrix analogue 
can be defined via the following Fourier representation: 

5 (Q - X^X) = f e-^2Trm-x^x)F^] dF„ = TT 17 (B.3) 

J 4tt 27r 

I Kj 

with the integration going over n x n real symmetric matrices: [Fn]ij = [Fn]ji- Substituting 
such a representation into the expression for Je{Q) and changing the order of integration 
over dFn and dX one may notice that the integral over X is essentially a product of A'^ 
identical Gaussian multivariable integrals ()A.6p where the role of A is played by the matrix 
el„ — iFn- The integrals are well-defined due to e > 0. Applying ()A.6p we arrive at 

Nn 

J,{Q) = {2-K) — Jn,N,e{Q)^ ^here 

and we have indicated explicitly the dependence on n and N for the sake of future reference. 
Notice that the integrand is invariant with respect to the rotations F — )• OFO~^ where O are 
orthogonal matrices satisfying 0^0 = 1. As Q is real symmetric matrix, it can be brought to 
the diagonal form by an orthogonal transformation. Hence the result of the integration can 
depend only on the eigenvalues qi,q2, ■■■,qn oi Q. Thus, it is enough to take Q to be diagonal 
from the very beginning. Now we separate the first eigenvalue from the rest: 

Q = diag(gi,g2,---,Q'n) = diag(q'i, Q„_i) 
and accordingly decompose the matrix Fn as 

where f = (/12, /13, fin) is a n — 1 component vector. 

Next step is to use the well-known property of the determinants composed of four blocks: 

det (el„ - iFn) = det (el„_i - iFn-i) (^e - ifn + f [el„_i - iFn-i]'^ f^^ 

which gives: 

Jn,NAQ) = j dFn-ie-i^^'i^"-'^"-') [det (el„_i - iFn-i)]-""^^ (B.6) 

' e - ^fn + f [eln-i - ^Fn^l] 
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The last integral over /n can be explicitly evaluated by using the formula 3.382.7 of 

where r(z^) is the Euler Gamma-function, and 0{x) = 1 for x > and zero otherwise. Taking 
into account e > 0, the result of the integration over /n gives 

^^^(.i) (I) exp {-^.1 (e + f - f^) } . (B.8) 

Now the integration over the vector df becomes the standard Gaussian and can be performed 
using (|A.6p yielding the factor: 



n-l 
1 \~ 



^ ' deti/2(el„_i-iF„_i^ 



2'Kqi 

Collecting all the factors we arrive at the recursive relation 



- JV-n-l 



Jn,NAQ) = 2nT[N/2) Ifj 9{q^)e--^^'^^Jn-i,N-iAQn-i) (B.9) 

This relation can be iterated further, and assuming N > n we arrive at the last step to (|B.7p 
which gives 

I N-n-l 

Jl,7V-n+l,.(gn)= ,^^^_,^,^ (f) ^ ^(g„)e-i^'^" (B.IO) 

and serves as an "initial condition" for our iteration scheme. This immediately yields the 
result: 

1 1 - - ^ 

2 2 TT 4 y^.^p r y-^) j=i 

Nn 

for N > n + 1. As Je{Q) = {2'r^)~ Jn,N ,e{Q) , in the limit e — )• the above relation yields 
precisely the required identity ()17|) . 



Appendix C. Parisi matrix, its eigenvalues and evaluation of traces in the 
replica limit. 

We start with describing the well known structure of the n x n matrix Q in the Parisi 
parametrisation, see Fig. 12a. At the beginning we set n diagonal entries Qaa all to the same 
value Qoia = 0. This value will be maintained at every but last step of the recursion. The 
off-diagonal part of the matrix Q in the Parisi scheme is built recursively as follows. At the 
first step we single out from the nxn matrix Q the chain consisting of n /mi blocks of the size 
mi < n, each situated on the main diagonal. All off-diagonal entries qap, a ^ j3 inside those 
blocks are filled in with the same value = Q'l < 0, whereas all the remaining n^(l — 1/mi) 
entries of the matrix Q are set to the value ^ < q^ < qi- The latter entries remain from 
now on intact to the end of the procedure, whereas some entries inside the diagonal nii x mi 
blocks will be subject to a further modification. At the next step of iteration in each of those 
diagonal blocks of the size mi we single out the chain of m^jmi smaller blocks of the size 
ni2 < nii^ each situated on the main diagonal. All off-diagonal entries (7^/3, a ^ j3 inside those 
sub-blocks are filled in with the same value qaj^ = q2 ^ Qi, whereas all the remaining entries of 
the matrix Q hold their old values. At the next step only some entries inside diagonal blocks 



Multifractality and Freezing 



31 



of the size m2 will be modified., etc. Iterating this procedure step by step one obtains after k 
steps a hierarchically built structure characterized by the sequence of integers 

n = nriQ > nii > m2 > • • • > ni^ > "ifc+i = 1 (CI) 

and the values placed in the diagonal blocks of the Q matrix satisfying: 

0<qo<qi<q2< ■■■<qk (C.2) 

Finally, we complete the procedure by filling in the n diagonal entries q^a of the matrix Q 
with one and the same value qaa = qd ^ qk- 

For the subsequent analysis we need the eigenvalues of the Parisi matrix Q. Those can 
be found easily together with the corresponding eigenvectors built according to a recursive 
procedure which uses the sequence Eq. (jC.ip . It is convenient to visualize eigenvectors as being 
"strings" of n boxes numbered from 1 to n, with l^^ component being a content of the box 
number I. 

At the first step i = 1 we choose the eigenvector to have all n boxes filled with the same 
content equal to unity. The corresponding eigenvalue is non-degenerate and equal to 

Ai = qd + qkimk-l) + qk-iimk-i-mk) + . . .+qi{mi-m2)+qoimo-mi){C.3) 

Now, at the subsequent steps i = 2, 3, . . . , A; + 2 one builds eigenvectors by the following 
procedure. The string of n boxes of an eigenvector belonging to i^^ family are subdivided 
into n/mj_i substrings of the length mj_i, and numbered accordingly by the index j = 
1, 2, . . . , n/mi_i. All mj_i boxes of the first substring j = 1 are filled invariably with all 
components equal to 1. Next we fill mj_i boxes in one (and only one) of the remaining 

— 1 substrings with all components equal to —1. In doing so we however impose a 

constraint that the substrings with the indices j given by j = 1 + ^ "^'"^ should be excluded 

from the procedure, with / being any integer satisfying 1 < / < — 1. After the choice 

of a particular substring is made, we fill all n — 2mj_i boxes of the remaining substrings 
with identically zero components. It is easy to see that all di = n/mi_i — n/mj_2 different 
eigenvectors of i^^ family built in such a way correspond to one and the same dj— degenerate 
eigenvalue 

Ai = qd+qkimk-l)+qk-i{mk-i-mk) + - ■ .+ft-i(?Tii_i -m^) -gi_2(mi_i)(C.4) 

In this way we find all n possible eigenvalues, the last being equal to 

Afc+2 = qd-qk mk+i = qd- qk- (C.5) 

The completeness of the procedure follows from the fact that sum of all the degeneracies di is 
equal to 

/ n \ f n n \ f n n \ 

1+ 1 + +•••+ 

\mi J \m2 niij \mk+i J 

Note that all the found eigenvalues are positive due to inequalities Eq. (IC.2p between various gj, 
which is required by the positive definiteness of the matrix Q. Note also that all eigenvectors 
built in this way are obviously linearly independent, although the eigenvectors belonging to 
the same family are not orthogonal. The latter fact however does not have any bearing for 
our considerations. 

To facilitate the subsequent treatment it is convenient to introduce the following 
(generalized) function of the variable g, see Fig. 12b: 

fe 

x{q) = n + y^(m;+i - mi) e{q - qi) (C.6) 
1=0 
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where we use the notation 9{z) for the Heaviside step function: 9{z) = 1 for z > and zero 
otherwise. In view of the inequahties Eq. (|C.l|C.2p the function x{q) is piecewise-constant 
non-increasing, and changes between n and 1 as follows: 

x{q < qo) =mo = n, x{qo < q < qi) = nii, . . . , x{qk^i < q < q^,) = ruk, x{q > qk) = nik+i = 1(C.7) 
Comparison of this form with Eq. ()C.6P makes evident the validity of a useful inversion formula: 



x{q) 



k 

- + E 

1=0 



1 



mi+i 



1 

nil 



Oiq - qi) 



(C.8) 



which will be exploited by us shortly. 

As observed by Crisanti and Sommers[37] one can represent the eigenvalues Eq. ()C.4p of 
the Parisi matrix in a compact form via the following remarkable identities: 



Ai 



Id 







Qd 



x{q) dq = nqo+ / x{q) dq, Ai+2 



qo 



Qd 



x{q)dq, i = 0,l,...,k{C.9) 



As a consequence, these relations imply for any analytic function g{x) the identity 

1 1 ^ 1 / f'^'^ \ /I 1 \ / Z''^'* \ 
-Tr[g{Q)] = -Tg{Xi)d, = -g(nqo+ x{q)dq)+T( )g( x{q)dq) (C.W) 

Next one observes that taking the derivative of the generalized function from Eq. (|C.8P produces 

k 



d_ 

dq 



1 



x{q) 



E 

1=0 



1 



mi+i 



— ) S{q-qi) 
mi/ 



This fact allows one to rewrite the sum in Eq. ()C.10P in terms of an integral, yielding 



-Tr [g{Q)] = -g ( nqo + 
n n 



Qd 



x{q)' 



Qo 



+ 



Qk+0 



90-0 



Qd 



xiq)dq] ^ 



(C.ll) 



dq, 



x{q) 



where the short-hand notation gibO designates the limit from below/above. Further performing 
integration by parts, and using x{q > q^) = 1, x{q < qo) = n, we finally arrive at 



^Tr[gm = - 
n n 



g\nqo + 



Qd 



x{q) dq 



Qo 



a 



Qd 



x{q) dq 



Qo 



+ 



Qk 



a 



QO 



Qd 



x{q) dq dq+g{qd-qk)-{C 



We are actually interested in the replica limit n — >• 0. According to the Parisi prescription in 
such a limit the inequality Eq. fjC.ip should be reversed: 



n 



< mi < m2 < . . . < m^ < ni^^i = 1 



(C.13) 



and the function x{q) is now transformed to a non-decreasing function of the variable q in the 
interval qo < q < qk, and satisfying outside that interval the following properties 

x{q < qo) = 0, and x{q > qk) = 1. (C.14) 

In general, such a function also depends on the increasing sequence of k parameters mi described 
in Eq.dMD • 

The form of Eq. (jaT2]l makes it easy to perform the limit n — )• explicitly, and to obtain 
after exploitation of Eq. (|62p an important identity Eq. (j64p helping to evaluate the traces in 
the replica limit. Finally, let us mention the existence of an efficient method of the "replica 
Fourier transform" allowing one to diagonalise (and otherwise work) with much more general 
types of hierarchical matrices, see |20] for more details. 
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^(ab),{ab) 



Appendix D. Stability analysis of the saddle-point solution 

Our starting point is the functional ^n{Q) from (jSOp whose extrema we look for in the 
space of positive definite matrices Q constrained to have the diagonal entries qaa = R^- 
The independent variables are all off-diagonal entries ^(Qf,) where {ab) stands for n(n — l)/2 
"ordered" pairs with a < b, and the stationary values are found from the equations ()A.3p . 
The stability matrix in this space is given by A(^ab),{cd) = 9g^~^^g^*^n(^) '^hich should be 
evaluated at the saddle-point solution. In a general situation we should distinguish three types 
of entries of that matrix: the diagonal entries 

the entries for the case when the ordered pairs {ab) and {cd) share one common replica, that 
is 

A^abUac) = [{Q^'),^ (Q-^),, + {Q-'),, iQ-'')J ,b<C (D.2) 

and a similar expression for ^(ab)^(cfe) , a < c, and finally the entries for the ordered pairs {ab) 
and {cd) which do not share any common replica: 

^(a.),M) = HQ-'),, iQ-'),, + iQ''),, {Q-')J , (D.3) 

If we are interested in investigating stability of the replica-symmetric solution, we should 

substitute to the above equations qab = 90)Va 7^ 6 as well as {Q~^),, = PdjVa and 

{Q~^),f^ = POjVa < b, with pd and po taken from ()55|56p . This gives for the entries of 
the stability matrix 

A(ab),iab) = M=pI+pI- 7 - go + 1)2 ' =M= POPd + pI (D.4) 

and ^(a6),{cd) = ^3 = 2pq. As discovered by De Almeida and Thouless [19j the 
eigenvalues/eigenvectors of such n{n — l)/2 x n(n — 1) /2 matrix can be found explicitly. There 
are three families of eigenvectors. The first family consists of a single "replica-symmetric" 
eigenvector ei with all components [eijj-^^^ = 1. The corresponding eigenvalue is equal to the 

sum of all entries in one row of A that is Ai = yli + 2(n — 2)^2 + ("'~^)^("~^) Next family 
consists of d2 = n — 1 eigenvectors c = — 1 with one replica index c singled out. For 

example, suppose that c = 1, then [e^] j-^^^ = if a = 1 or 6 = 1, and [eg] ^-^^^ = — 1 otherwise 
(note that such eigenvector is orthogonal to ei). The corresponding eigenvalue shared by all 
the eigenvectors in the family is X2 = Ai + {n — 4)^2 — (^^ — 3)^3 . Finally , third family consists 
of da = "'*'"g eigenvectors e^^'^^ with an ordered pair of replica indices c < d singled out. For 



example, if {cd) = (12) then components of the corresponding eigenvector are 



(12) 



,(11) 



(12) 



= if a = 1, 2 or 6 = 1, 2, and otherwise 

(ab) 



(11) 



= p where the values of ^, V', p 

(ab) 



should be chosen to make e^^^^ orthogonal to e'^^ and ei. The eigenvalue shared by the third 
family turns out to be given by an n— independent expression A3 = — 2^2 + ^3- Since 
1 + (^2 + da = n{n — l)/2 no more eigenvalues are possible. 

It is well known in general (and can be easily checked for our model) that it is third 
family which gives rise to "dangerous" fluctuations breaking down the replica symmetry of 
the saddle-point solution in the limit n — )• below some critical temperature at which A3 
vanishes. Substituting the expressions ()D.4p into A3 and using the relation ([57|) we find that 
the condition A3 = is equivalent to 
1 7I/2 



R^-qo R^-qo + 1 



(D.5) 
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Finally, using for the combination R'^ — qq = cIq the equation ([^^ we find after simple algebra 
that the critical value of the parameter = f3g is given by ^/7^ = as was quoted in the 
text. 

Let us now turn to the stability issue for the one-step RSB solution which we claimed 
to be the correct choice below the critical temperature. According to Fig. 12a the one-step 
solution is characterized by the matrices Q with all diagonal entries still equal to R^, and two 
different values of the off-diagonal entries qi > qo- The size of blocks containing qi is equal to 
m. Such more complicated structure of Q generates more types of different elements in the 
stability matrix A, and although its eigenvalues/eigenvectors can be still successfully found 
[16], actual analysis becomes long. Referring the interested reader to Appendix B3 of [16] for a 
detailed exposition, we give below a very brief summary of the the outcome of the procedure. 
Actually the stability analysis of |16j was performed for models with general random potential 
characterized by the covariance (xi) F (X2)) = Nf ((xi — X2)^/2iV). It was found that 
the matrix A in that general case has nine different eigenvector families, and the stability is 
controlled by two of them with eigenvalues given by 



A^ 



and 



1 



1 



R"^ -qi + m{qi - qo) 



Vf"iR'-Qo) 



1 



1 



R'^ -qi+ m{qi - qo] 



qo) 



(D.6) 



A^ 



K 



R^ 



qi 



qij 



qi 



qi) 



(D.7) 



assuming f"{x) > 0. If both of the above eigenvalues are non-negative, and the second 
derivative f"{x) is monotonically decreasing with x then all the remaining seven eigenvalues 
are strictly positive and the system is stable. 

It is easy to understand that the logarithmic case (llip considered in this paper is 



recovered, after all due rescalings in the limit ^ 1, by the choice f{x) = — (7 ln(x + l) 
so that f"{x) = Using now the equilibrium values for the parameters qi,qQ,m given 
in (|77p one finds after a straightforward algebra that the first brackets in ()D.6p . ()D.7p identically 
vanish leaving us with Aq = A^ = everywhere in the low-temperature phase. This implies 
indeed that the corresponding one-step RSB solution is marginally stable. 
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